An EOF Iteration Approach for Obtaining Homogeneous Radiative 
Fluxes from Satellites Observations 

Banglin Zhang (1), Rachel T. Pinker (2) and Paul W. Stackhouse, Jr. (3) 

(1) Global Modeling and Assimilation Office, Code 610.1 
SAIC, Goddard Space Flight Center Greenbelt, MD 20771 
(2) Department of Meteorology 
University of Maryland, College Park, MD 20742 
(2) NASA Langley Research Center 
2 1 Langley Boulevard, Mail Stop 420 
Hampton, VA 23681-2199 
Submitted to Journal of Applied Meteorology 
September 4, 2005 
Revised version: 

May 2006 
Final Version 
August 2006 
Corresponding Author: 

Rachel T. Pinker 
University of Maryland 

Department of Atmospheric and Oceanic Science 
College Park, MD 20742 

pinker@atmos.umd.edu Tel: 301-405-5380 FAX: 301-314-9482 


1 



Abstract 


Conventional observations of climate parameters are sparse in space and/or in 
time and the representativeness of such infonnation needs to be optimized. Observations 
from satellites provide improved spatial coverage than point observations however they 
pose new challenges for obtaining homogeneous coverage. Surface radiative fluxes, the 
forcing functions of the hydrologic cycle and biogeophysical processes, are now 
becoming available from global scale satellite observations. They are derived from 
independent satellite platforms and sensors that differ in temporal and spatial resolution 
and in the size of the footprint from which information is derived. Data gaps, degraded 
spatial resolution near boundaries of geostationary satellites, and different viewing 
geometries in areas of satellite overlap, could result in biased estimates of radiative 
fluxes. In this study, discussed will be issues related to the sources of inhomogeneity in 
surface radiative fluxes as derived from satellites; development of an approach to obtain 
homogeneous data sets; and application of the methodology to the widely used 
International Satellite Cloud Climatology Project (ISCCP) data that currently serve as a 
source of infonnation for deriving estimates of surface and top of the atmosphere 
radiative fluxes. Introduced is an Empirical Orthogonal Function (EOF) iteration scheme 
for homogenizing the fluxes. The scheme is evaluated in several ways including 
comparison of the inferred radiative fluxes against ground observations, both before and 
after the EOF approach is applied. On the average, the latter reduces the rms error by 
about 2-3 W/m 2 . 
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1 . 


Introduction 


1.1 At issue 

To advance the understanding of the water cycle and land-atmosphere interactions 
on global scale, information on radiative fluxes is needed at similar scales, and in principle, 
can be obtained from satellite observations. Only recently global scale satellite 
observations at climatic time scales have become available. While issues related to paucity 
of data have been addressed in the past, each climatic parameter poses a unique challenge 
for obtaining homogeneous time series representative of large spatial scales. The most 
widely used parameter to study climate variability is surface temperature, either at shelter 
level over land stations, or sea surface temperature (SST), as derived both from satellites 
and in-situ observations. Reynolds and Smith (1995) developed an optimal interpolation 
(01) technique implemented weekly at 1° resolution using ship and buoy data of SST and 
satellite data. In their approach attention is given to the removal of large-scale biases of 
satellite data, as compared to in situ data. Smith et al. (1996) expanded this approach by 
using the best estimates of the full covariance structure of the mean SST fields to fit the in 
situ data and fill data gaps. They used spatial functions defined by empirical orthogonal 
functions (EOF), computed from a principal component analysis (PCA). It is claimed that 
this approach is an improvement over the 01 method, which utilized local interpolation 
with modeled covariance functions, while in the newer scheme basin wide covariance 
structures are used. Other studies dealt with various aspects of interpolation needs of 
climatic parameters, emphasizing different aspects of the problem. For example, Rayner et 
al. (1995) attempt to improve global long tenn sea ice and SST analysis while Shriver and 
O’Brien (1995) address the wind stress issue. Radiative fluxes as derived from satellites 
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are more recent climate parameters that are becoming available at spatial and temporal 
scales of interest in climate research. These pose a different type of difficulties than those 
addressed previously, and will be addressed in this study. 

1.2 Satellite data to be used 

The International Satellite Cloud Climatology Project is a major program of the 
World Climate Research Program (WCRP) for preparing satellite observations, which are 
suitable for inferring global distribution of clouds, their properties, and their effect on the 
radiation balance (Schiffer and Rossow, 1983; Schiffer and Rossow, 1985; Rossow and 
Schiffer, 1991; Rossow and Garder, 1993a and 1993b; Brest et al., 1997; Rossow and 
Schiffer, 1999). Data collection began on July 1, 1983, continued through 2006 and there 
are plans for extension beyond 2006. 

The ISCCP D1 data are at 2.5 degree spatial resolution and 3-hourly time intervals. 
They are based on observations made from several geostationary and polar orbiting 
satellites. These data have been extensively used for inferring global scale distributions of 
clouds, cloud properties, and surface radiative fluxes (Han et al. 1995; Zhang et al. 1995; 
Pinker and Laszlo, 1992; Tselioudis et al. 2000; Zhang et al. 2004; Pinker et al. 2005). 
Geophysical parameters derived from the ISCCP D1 data have also been used in a wide 
range of applications (Hahmann et ah, 1995; Bony et al. 1997; Machado et al. 1998; Prigent 
et al. 1998; Sui et al. 2003). For studies of the hydrologic cycle at regional scale, 
information on surface radiative fluxes is needed at high spatial resolution. The ISCCP DX 
satellite observations (the basis for the D1 data-set) are made at 30 km resolution from 
several geostationary and polar orbiting satellites as provided by various agencies. The 
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Tiros-N data come from the National Oceanic and Atmospheric Administration (NOAA); 
the METEOSAT data are provided by the European Meteorological Satellite Organization 
(EUMETSAT); the GOES-EAST data are provided by the Meteorological Service of 
Canada; the GOES-WEST data come from the University of Colorado, Cooperative 
Institute for Research in the Atmosphere (CIRA); and the GMS data come from the Japan 
Meteorological Agency. The observations are not homogeneous in terms of spatial and 
temporal coverage; in particular, when geostationary satellites are used, merging issues 
arise. 


1.2.1 Motivation for study 

Data from geostationary satellites at the ISCCP DX level are provided 
independently for each satellite. An ISCCP DX product is also available from polar 
orbiting satellites, and is provided on an orbit-by-orbit basis, covering the entire globe. 
When all available observations at the DX level are combined from all the observing 
platforms, the coverage density is as presented in Figure 1 (a). At any particular time 
and location, climatic information inferred from these observations will be affected by 
the uneven distribution of the observations. The objective of this study is to develop an 
approach for optimal smoothing and interpolation of the derived surface radiative fluxes 
as obtained from the raw satellite observations. Discussion will focus on the surface 
short-wave radiative fluxes derived from the ISCCP DX data from GOES and 
METEOSAT satellites only, since in these regions ample ground observations are 
available for the evaluation of the methodology. 
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1.2.2 Methodology to derive surface radiative fluxes 


The primary objective of this study is to develop a methodology that will provide 
homogeneous infonnation on surface and top of the atmosphere radiative fluxes for large 
areas. For estimating the surface short-wave (SW), Photosynthetically Active Radiation 
(PAR), and near-infrared (NIR) fluxes from the satellite observations at both boundaries 
of the atmosphere from the ISCCP DX data, an inference technique developed at the 
University of Maryland (Pinker and Laszlo 1992) is used. The method is based on 
radiative transfer theory, and produces direct and diffuse fluxes in five spectral intervals 
in the 0.2-4. 0 pm range. The radiative transfer model accounts for the absorption and 
scattering processes occurring in the atmosphere and for the interaction of the radiation 
with the surface. The interaction with the surface is modeled spectrally, and allows 
incorporation of spectral observations that are becoming available from EOS 
observations (MODerate resolution Imaging Spectrometer (MODIS) (King et ah, 1992) 
and the Global Imager (GLI) sensor, on the Advanced Earth Observing Satellite 
(ADEOS-II)). 

The radiative fluxes at the boundaries of the atmosphere are computed by 
determining atmospheric transmission and reflection (optical functions) and the surface 
albedo, pertaining to a particular satellite observation. The retrieved optical functions, 
along with the surface albedos, are then used to independently compute the fluxes for 
clear and cloudy conditions. The all-sky flux is obtained by using infonnation on cloud 
cover. The fluxes for this study were produced at the University of Maryland with 
Version 2.1 of the Global Energy and Water cycle EXperiment/Surface Radiation Budget 
(UMD/SRB) model (Whitlock et al. 1995; Pinker et al. 1995), which includes improved 


6 



parameterization of water vapor absorption (Ramaswamy and Freidenreich, 1992) as well 
as improved techniques for replacing missing observations as compared to Version 1.1. 
Missing values were replaced based on time-space interpolation. This latter version has 
been implemented at the Surface Radiation Budget/Satellite Data Analysis center 
(SDAC) at NASA Langley Research Center with ISCCP Cl data at 2.5 degree resolution 
for July 1983-December 1988, and results are available from the Langley DAAC at: 
http://eosweb.larc.nasa.gov/ . More recently, it is being implemented at LaRC at 1 degree 
resolution and results are available at: http://srb- 

swlw.larc.nasa.gov/GEWEX_SRB_homepage.html 

The monthly mean shortwave downward flux (W/m“) for January 1992 as derived 
from ISCCP DX data from GOES, METEOSAT, and GMS is illustrated in Figure 2 (b). 
To obtain homogeneous and coherent radiative fluxes from the satellite observations 
requires addressing issues related to uneven sampling, data gaps, degraded spatial 
resolution as one moves away from the nadir, and different viewing geometries in areas 
of overlap. The Empirical Orthogonal Function (EOF) iteration scheme is well suited for 
getting integrated, smoothed, high-resolution surface radiative fluxes, because in this 
approach, spatial modes are derived statistically from available observations and then 
used to reconstruct the information in areas of degraded data quality. It should be pointed 
out that the ISCCP data are normalized to one reference polar orbiting satellite 
(afternoon) measurements and in addition, the residuals from the normalization are 
further corrected by calibration adjustment, so that certain inhomogeneity/incoherence in 
the ISCCP DX data are removed to begin with. Selected available interpolation 
methodologies and the one used in this study will be presented in Section 2. Evaluation 
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of the proposed methodology will be discussed in Section 3. Conclusions and discussion 
will be presented in Section 4. 

2. Interpolation schemes 

2.1 Available methodologies 

For data gap filling and interpolation, several schemes have been developed. For 
example, Renka (1997a, 1997b) developed an interpolation approach using piecewise 
linear interpolation functions. Optimal interpolation (01) techniques are more widely 
used for meteorological applications (Lorenc, 1981). 01 is performed locally with 

modeled covariance functions that take into account local data structure. The 
interpolated fields from both 01 and simple interpolation schemes are usually used for 
relatively small regions and can result in isolated anomalies due to unreasonable 
stretching of data to fill the gaps. 

Smith et al. (1996) developed an interpolation method using EOFs from 01 
analyses to remove deficiencies in SST interpolated fields. The core of the method is 
based on the adaptation of EOF reconstruction techniques used to fill gaps in the gridded 
data-sets. The approach involves (i) development of spatial EOFs from recent data-rich 
years; (ii) use of dominant EOF modes as basis functions that are fitted to available data 
from earlier data-poor years; and (iii) reconstruction of global SST fields from these 
spatial and temporal modes. The EOFs efficiently use the best estimate of the full 
covariance structure of the monthly mean SST fields to fit in situ SST observations to 
data voids and field smoothing. This method can be used to interpolate other 
atmospheric and oceanograghic data. 
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2.2 Empirical Orthogonal Function (EOF) iteration scheme 
When applied to satellite estimates of radiative fluxes that have both missing and 
degraded data, the Smith et al. EOF reconstruction scheme can not distinguish between 
good and degraded quality observations. This problem can be addressed by using an 
EOF iteration approach. The scheme was first developed by Zhang et al (1993) to 
forecast multi-dimensional time series and was applied to ENSO prediction (Zhang and 
Pan, 1996). Before applying the EOF iteration scheme for interpolation and gap filling, 
the data are first divided into two groups: a) high quality, such as observations near nadir 
where the frequency of observations is high (anchor region); b) lower quality, such as the 
regions near viewing edges of each satellite; regions of overlap with different viewing 
geometry; and regions where data are missing. In Figure 2 illustrated are two examples 
of possible selection of anchor regions for ISCCP DX GOES and METOSAT satellites. 
In order to derive a “first guess” field for the EOF analysis, used is Version 2.1 of 
UMD/SRB model to process at daily time scale the ISCCP DX data from METEOSAT 
and GOES satellite for the period 1989-1993 at 0.5-degree resolution to derive the daily 
radiative fluxes. For each month, the daily fluxes in areas of overlap are averaged and 
missing data are replaced with fluxes derived from the ISCCP D1 data, interpolated to 
0.5 -degree spatial resolution, using the interpolation approach developed by Renka 
(1997a, 1997b). These data serve as a first guess in the iteration process. We create a 
data matrix, which consists of the departures of the “first guess” field from the monthly 
mean, in the form: 
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where N is the number of time series (for example, it is 3 1 daily average values of surface 
short-wave flux deviations from the January monthly mean), and M is the number of all 
spatial grids for the data set. 

The EOF analysis is applied on the data matrix F (0) to produce eigenvalues in a 
descending order (from larger to smaller ones) and associated EOF modes. The first EOF 
accounts for most of the total variance. Each following EOF explains a smaller part of 
the total variance. The ratio between the accumulated eigenvalues of the first K EOFs 
and the sum of all eigenvalues represents the percentage variance explained and is often 
the basis for deciding what should be the number of EOFs to retain. For example, in 
Figure 3, illustrated is the explained variance as a function of the number of EOF modes 
used for January and July of 1992. We use the ratio of 80% to select the number of EOF 
modes which explain 80% of the total variance of the first guess field. In this case, 
twenty EOFs will explain about 80 % of the total variance. Therefore, the first twenty 
EOFs are selected to reconstruct the data in the first iteration step. Subsequently, the 
EOF analysis is repeated using the available data in the high quality regions (the “anchor” 
regions), and the reconstructed data in the poor quality region. The resulting “fusing” 
scheme will iteratively replace the data in the degraded domain by propagating 
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information from the high quality regions. It should be pointed out that the EOF iteration 
scheme does not address problems that may be caused by biases between different 
instruments. 

The initial reconstructed field is obtained as follows: 
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Here yT and xT are the principal components and EOF modes derived from the first 
guess Fnm • Subsequently, the missing values and degraded data are replaced with the 
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This reconstructed field is then used to replace the missing and degraded values to create 
F ( nxm for next EOF analysis. Iteratively, the EOF analysis is repeated up to the n th step, 
the degraded data and missing values are replaced each time until the process converges. 
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In the final step, EOF smoothing is applied on the entire data set, including the “anchor” 
regions, resulting in the final reconstructed and smoothed field given as: 
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In Figure 4, the convergence process is illustrated for January and July of 1992. The 
RMS difference for the entire field between successive iteration steps is calculated, 
nonnalized by the standard deviation. When this difference is less than 0.5 %, the 
iteration process is complete. Figure 4 shows that about 7-8 iterations are needed for 
converge. The final filled values are the sum of the base monthly mean and the 
interpolated values from the EOF iteration scheme. 

3. Evaluation of the merging process 
3.1 Qualitative evaluation 

The EOF iteration scheme was applied to surface shortwave radiation derived 
from METEOSAT and GOES at daily time scale for January 1992, over (70° S -70° N, 
180° W-70° E). Illustrated are the results for the two control tests. In the first test, the 
sub-region 20° S-20° N, 130°-70° W where observations are most dense is chosen as a 
test region, and the computed fluxes are removed for the entire month from all grids 
except for 5 by 5 grids. For the artificially missing values, the Renka scheme is used to 
get a first guess, and subsequently, the EOF scheme is run with and without the iteration 
approach (without the iteration it is equivalent to the Smith method). The daily 
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shortwave downward fluxes on January 2, 1992 are shown in Figure 5 (a) for actual 
observations, as estimated by the Renka interpolation scheme, the Smith EOF 
reconstruction scheme, and EOF iteration scheme. No significant difference is observed 
for the three schemes. The RMS differences between actual daily values and estimations 
are 30.13, 27.19, and 27.43 Wm" 2 for Renka scheme, Smith scheme, and EOF iteration 
scheme, respectively. 

In the second control test, data are removed from the entire sub-region 20 S-20 N, 
130-70 W on even days for the month of January 1992 Figure 5 (b) shows the original 
field before removing data for the tests and estimates from the three interpolation 
schemes on January 2, 1992. As seen, the results from all 3 interpolation schemes have 
degraded when compared with the control test. However, results from the EOF iteration 
scheme are smoother than the results from the Renka or Smith schemes. The RMS 
differences between the original daily values and estimates are 53.66, 49.74, and 47.36 
Wm' for the Renka, Smith, and EOF iteration schemes, respectively. Based on both 
spatial structure and RMS differences, the EOF iteration scheme seems to perform better 
than the Renka and Smith schemes. 

Sensitivity tests were conducted by using two possible ways to select “anchor” 
regions and regions of degraded information as illustrated in Figure 2 for January 1992. 
It was shown that the differences in the derived fluxes using the different anchor regions 
were quite small. This was repeated for July 1992 and again, the differences in the 
derived fluxes were small. 
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3.1 Quanitative evaluation 

A quantitative evaluation of the impact of the proposed merging methodology on 
the derived surface fluxes was performed by comparison of results from Renka scheme, 
Smith method, and EOF iteration scheme against ground observations on a daily time 
scale. For this purpose, used was a partial archive of ground truth as available at the 
NASA Fangley Research Center and as described in Chiacchio et al. (2004). The 
broadband SW surface measurements were provided by the World Radiation Data Center 
(WRDC/GEBA)/Canadian Networks. Chosen were sites with continuous long-term 
observations in regions observed by the two geostationary satellites GOES and 
METEOSAT. Used was the period 1989 to 1993 for which the proposed methodology 
has been implemented. The sites used in this study are shown in Figure 6. The analysis 
was done separately for North America and for Africa. Used were all available stations 
for each region which resulted in 9447 data points for North America and 59796 data 
points for Africa; independent statistics were performed for several sub-regions clustered 
by proximity to each other, as show in Figure 8. The complete statistics for each sub- 
region as well as for North America and for Africa is presented in Table 1 . As evident, 
when comparing Renka method, the Smith method and EOF iteration scheme reduced the 
average rms was reduced by about 2-3 Wm' while the bias was reduced by about 5-6 
Wm' . The bias and rms averages from EOF iteration are slightly smaller than from the 
Smith method. It is believed that in regions of highly degraded observations or sparse 
coverage, the improvements would be better. Due to the sparcity of observations in such 
regions it was not possible to undertake such specific analysis at this time. Illustration of 
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results for all N. America using all ground observations and for Africa, using 
observations in region III only, is given in Figure 7. 

4, Discussion 

The inhomogeneity of satellite observations used in the evaluation of estimated 
surface parameters (e.g., surface radiation or surface temperature) against ground 
observations will impact the demonstrated accuracy. Similar problems arise when 
satellite estimates are used to evaluate climate or numerical weather prediction models, 
which have homogeneous grids. Therefore, it is important to ‘homogenize’ the various 
satellite-based parameters in a way that the effect of inhomogeneity and incoherence will 
be minimized. In the present study, an effort in that direction has been made. 
Specifically, an EOF analysis was used to utilize observations in regions of frequent 
observations, to propagate temporal and spatial structural information into regions of 
missing or degraded observations. It was demonstrated that this approach resulted in 
improved estimates of surface fluxes as well as in smoothed and coherent fields. 
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The frequency distribution of ISCCP geostationary and polar orbiting 
satellite observations at DX resolution for October 1986. 

Monthly mean shortwave downward flux (W/m ) for January 1992 
derived from ISCCP DX GOES, METEOSAT, and GMS observations. 

Illustration of two (top and bottom) possible ways to select anchor regions 
(rectangle in the middle) for EOF iteration. Inside rectangle: anchor 
region; outside rectangle: degraded region; hatched: overlap region; white: 
missing values. Upper: larger anchor region; Lower: smaller anchor 
region. 

The explained % variance of shortwave surface fluxes by EOFs, as a 
function of the number of EOF modes used for January (top) and July 
(bottom) of 1992. 

The RMS difference, normalized by the total variance, between successive 
iteration steps. When this difference becomes less than 0.5 %, the 
iteration process is stopped. 

The daily shortwave downward fluxes for control test 1 on January 2, 
1992 from: original observations (upper left panel); estimated by Renka 
interpolation scheme (lower left panel); Smith’s EOF reconstruction 
scheme (upper right panel); and proposed EOF iteration scheme (lower 
right panel). 

Same as Figure 6 but for control test 2. 
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Figure 6. Location of stations used for evaluation of the improvements in surface 
radiative fluxes as a result of the proposed merging procedure. Circle: 
North America Region 1 (NA1); oentagram: North America Region 2 
(NA2); hexagram: Africa Region 1 (AF1); asterick: Africa Region 2 (AF2); 
square: Africa Region 3 (AF3); and diamond: Africa Region 4 (AF4). 

Figure 7. Comparison of merged surface SW radiative fluxes as derived from 
ISCCP DX data gridded to 0.5° against ground observations for the period 
1989-1993. Over North America, results are illustrated for all available 
stations. Over Africa, results are illustrated only for stations in region III. 


22 




BON 
60N 
4-0 N 
20 N 
EQ 
20S 
40S 
60S 
80S 


40 GO BO lOO 


120 ISO ISO 24-0 300 360 420 480 340 



Figure 1. Upper: The frequency distribution of ISCCP geostationary and polar 

orbiting satellite observations at DX resolution for October 1986. Lower: 
Monthly mean shortwave downward flux (W/m ) for January 1992 
derived from GOES, METEOSAT, and GMS ISCCP DX data 
independently. 


23 



60 N 
40 N 
20N 
0 

20S 
40S 
60S 

180 1 20W 60W 0 60E 

Figure 2. Illustration of two (top and bottom) possible ways to select anchor regions 
(rectangle in the middle) for EOF iteration. Inside rectangle: anchor 
region; outside rectangle: degraded region; hatched: overlap region; white: 
missing values. Upper: larger anchor region; Lower: smaller anchor 
region. 
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Figure 5. Upper four panels: The daily shortwave downward fluxes over region for control 
test 1 on January 2, 1992 from actual observations (upper left panel), estimated 
by Re nk a interpolation scheme (lower left panel), Smith’s EOF reconstruction 
scheme (upper right panel), and EOF iteration scheme (lower right panel). 
Lower four panels: Same as upper but for control tests 2. 
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Station Locations 



Figure 6. Location of stations used for evaluation of the improvements in surface 
radiative fluxes as a result of the proposed merging procedure. Circle: 
North America Region 1 (NA1); oentagram: North America Region 2 
(NA2); hexagram: Africa Region 1 (AF1); asterick: Africa Region 2 (AF2); 
square: Africa Region 3 (AF3); and diamond: Africa Region 4 (AF4). 
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Figure 7. Comparison of merged surface SW radiative fluxes as derived from 
ISCCP DX data gridded to 0.5° against ground observations for the period 
1989-1993. Over North America, results are illustrated for all available 
stations. Over Africa, results are illustrated only for stations in region III. 
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Table 1. Evaluation of the Ranka, Smith and the EOF merging scheme using ground 
observations of radiative fluxes. Given are the bias numbers (first) and the 
rms (second). When all observations over N. America and Africa were 
used, two experiments were performed: a) using all data; b) using data from 
areas where no satellite observations were available and the satellite values 
used are the result of the EOF fusing method. 


Regions 

Ranka Interpolation 

Smith Method 

EOF M 

lerging 


All data 

Missing 

Data 

All data 

Missing 

Data 

All data 

Missing 

Data 

N America 

-8.41/38.83 

-12.06/52.66 

-8.08/38.84 

-6.62/52.52 

-7.50/36.06 

-6.61/48.83 

Africa 

-13.29/44.83 

-45.07/69.68 

-10.63/41.73 

-43.53/61.73 

-10.65/40.96 

-43.50/61.62 

N America I 

-5.32/37.67 


-5.62/36.99 


-4.88/35.56 


N America II 

-11.56/37.41 


-10.90/37.50 


-10.16/33.56 


Africa I 

-22.35/43.18 


-20.42/41.17 


-20.42/40.27 


Africa II 

-1.77/42.81 


-0.99/40.91 


-0.99/39.71 


Africa III 

-7.29/33.07 


-4.64/30.15 


-4.67/29.74 


Africa IV 

-19.57/45.97 


-19.48/43.57 


-19.53/43.29 
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